a. Field of the Invention
The invention relates generally to computerized techniques for processing data obtained from radar to track multiple discrete objects.
b. Description of the Background
There are many situations where the courses of multiple objects in a region must be tracked. Typically, radar is used to scan the region and generate discrete images or xe2x80x9csnapshotsxe2x80x9d based on sets of returns or observations. In some types of tracking systems, all the returns from any one object are represented in an image as a single point unrelated to the shape or size of the objects. xe2x80x9cTrackingxe2x80x9d is the process of identifying a sequence of points from a respective sequence of the images that represents the motion of an object. The tracking problem is difficult when there are multiple closely spaced objects because the objects can change speed and direction rapidly and move into and out of the line of sight for other objects. The problem is exacerbated because each set of returns may result from noise as well as echoes from the actual objects. The returns resulting from the noise are also called false positives. Likewise, the radar will not detect all echoes from the actual objects and this phenomena is called a false negative or xe2x80x9cmissed detectxe2x80x9d error. For tracking airborne objects, a large distance between the radar and the objects diminishes the signal to noise ratio so the number of false positives and false negatives can be high. For robotic applications, the power of the radar is low and as a result, the signal to noise ratio can also be low and the number of false positives and false negatives high.
In view of the proximity of the objects to one another, varied motion of the objects and false positives and false negatives, multiple sequential images should be analyzed collectively to obtain enough information to properly assign the points to the proper tracks. Naturally, the larger the number of images that are analyzed, the greater the amount of information that must be processed.
While identifying the track of an object, a kinematic model may be generated describing the location, velocity and acceleration of the object. Such a model provides the means by which the future motion of the object can be predicted. Based upon such a prediction, appropriate action may be initiated. For example, in a military application there is a need to track multiple enemy aircraft or missiles in a region to predict their objective, plan responses and intercept them. Alternatively, in a commercial air traffic control application there is a need to track multiple commercial aircraft around an airport to predict their future courses and avoid collision. Further, these and other applications, such as robotic applications, may use radar, sonar, infrared or other object detecting radiation bandwidths for tracking objects. In particular, in robotic applications reflected radiation can be used to track a single object which moves relative to the robot (or vice versa) so the robot can work on the object.
Consider the very simple example of two objects being tracked and no false positives or false negatives. The radar, after scanning at time t1, reports objects at two locations in a first observation set. That is, the radar returns a set of two observations {o11, o12}. At time t2 the radar returns a similar set of two observations {o21, o22) from a second observation set. Suppose from prior processing that track data for two tracks T1 and T2 includes the locations at t0 of two objects. Track T1 may be extended through the points in the two sets of observations in any of four ways, as may track T2. The possible extensions of T1 can be described as: {T1, o11, o21}, {T1, o11, o22}, {T1, o12, o21} and {T1, o12, o22}. Tracks can likewise be extended from T2 in four possible ways, including {T2, o12, o21 }. FIG. 1 illustrates these five (out of eight) possible tracks (to simplify the problem for purposes of explanation). The five track extensions are labeled h11, h12, h13, h14, and h21wherein h11 is derived from {T1, o11, o21}, h12 is derived from {T1, o11, o22}, h13 is derived from {T1, o12, o21}, h14 is derived from {T1, o12, o22}, and h21 is derived from {T2, o11, o21}. The problem of determining which such track extensions are the most likely or optimal is hereinafter known as the assignment problem.
It is known from prior art to determine a figure of merit or cost for assigning each of the points in the images to a track. The figure of merit or cost is based on the likelihood that the point is actually part of the track. For example, the figure of merit or cost may be based on the distance from the point to an extrapolation of the track. FIG. 1 illustrates costs xcex421 and xcex422 for hypothetical extension h21 and modeled target characteristics. The function to calculate the cost will normally incorporate detailed characteristics of the sensor, such as probability of measurement error, and track characteristics, such as likelihood of track maneuver.
FIG. 2 illustrates a two by two by two matrix, c, that contains the costs for each point in relation to each possible track. The cost matrix is indexed along one axis by the track number, along another axis by the image number and along the third axis by a point number. Thus, each position in the cost matrix lists the cost for a unique combination of points and a track, one point from each image. FIG. 2 also illustrates a {0, 1} assignment matrix, z, which is defined with the same dimensions as the cost matrix. Setting a position in the assignment matrix to xe2x80x9conexe2x80x9d means that the equivalent position in the cost matrix is selected into the solution. The illustrated solution matrix selects the {h14, h21} solution previously described. Note that for the above example of two tracks and two snapshots, the resulting cost and assignment matrices are three-dimensional. As used in this patent application, the term xe2x80x9cdimensionxe2x80x9d means the number of axes in the cost or assignment matrix while size refers to the number of elements along a typical axis. The costs and assignments have been grouped in matrices to facilitate computation.
A solution to the assignment problem satisfies two constraintsxe2x80x94first, the sum of the associated costs for assigning points to a track extension is minimized and, second, if no false positives or false negatives exist, then each point is assigned to one and only one track.
When false positives exist, however, additional hypothetical track extensions incorporating the false positives will be generated. Further note that the random locations of false positives will, in general, not fit well with true data and such additional hypothetical track extensions will result in higher costs. Also note that when false negative errors exist, then the size of the cost matrix must grow to include hypothetical track extensions formulated with xe2x80x9cgapsxe2x80x9d (i.e., data omissions where there should be legitimate observation data) for the false negatives. Thus, the second criteria must be weakened to reflect false positives not being as signed and also to permit the gap filler to be multiply assigned. With hypothetical cost calculated in this manner then the foregoing criteria for minimization will to materialize the false negatives and avoid the false positives.
For a three-dimensional problem, as is illustrated in FIG. 1, but with N1 (initial) tracks, N2 observations in scan 1, N3 observations in scan 2, false positives and negatives assumed, the assignment problem can be formulated as:                                                         (              a              )                                                          Minimize:                                                                          ∑                                                      i                    1                                    =                  0                                                  N                  1                                            ⁢                                                ∑                                                            i                      2                                        =                    0                                                        N                    2                                                  ⁢                                                      ∑                                                                  i                        3                                            =                      0                                                              N                      3                                                        ⁢                                                            C                                                                        i                          1                                                ⁢                                                  i                          2                                                ⁢                                                  i                          3                                                                                      ⁢                                          z                                                                        i                          1                                                ⁢                                                  i                          2                                                ⁢                                                  i                          3                                                                                                                                                                            xe2x80x83                                                                          (              b              )                                                          Subject   to:                                                                                                                ∑                                                                  i                        2                                            =                      1                                                              N                      2                                                        ⁢                                                            ∑                                                                        i                          3                                                =                        1                                                                    N                        3                                                              ⁢                                          z                                                                        i                          1                                                ⁢                                                  i                          2                                                ⁢                                                  i                          3                                                                                                                    =                1                            ,                                                                                            i                  1                                =                1                            ,                              …                ⁢                                  xe2x80x83                                ⁢                                  N                  1                                            ,                                                                          (              c              )                                                          xe2x80x83                                                                                                                ∑                                                                  i                        1                                            =                      1                                                              N                      1                                                        ⁢                                                            ∑                                                                        i                          3                                                =                        1                                                                    N                        3                                                              ⁢                                          z                                                                        i                          1                                                ⁢                                                  i                          2                                                ⁢                                                  i                          3                                                                                                                    ≤                1                            ,                                                                                            i                  2                                =                1                            ,                              …                ⁢                                  xe2x80x83                                ⁢                                  N                  2                                            ,                                                                          (              d              )                                                          xe2x80x83                                                                                                                ∑                                                                  i                        1                                            =                      1                                                              N                      1                                                        ⁢                                                            ∑                                                                        i                          2                                                =                        1                                                                    N                        2                                                              ⁢                                          z                                                                        i                          1                                                ⁢                                                  i                          2                                                ⁢                                                  i                          3                                                                                                                    ≤                1                            ,                                                                                            i                  3                                =                1                            ,                              …                ⁢                                  xe2x80x83                                ⁢                                  N                  3                                            ,                                                                          (              e              )                                                          xe2x80x83                                                                          z                                                      i                    1                                    ⁢                                      i                    2                                    ⁢                                      i                    3                                                              ∈                              {                                  0                  ,                  1                                }                                                                                        ∀                                  xe2x80x83                                ⁢                                  z                                                            i                      1                                        ⁢                                          i                      2                                        ⁢                                          i                      3                                                                                  ,                                                          (        0.1        )            
where xe2x80x9ccxe2x80x9d is the cost and xe2x80x9czxe2x80x9d is a point or observation assignment, as in FIG. 2.
The minimization equation or equivalently objective function (0.1)(a) specifies the sum of the element by element product of the c and z matrices. The summation includes hypothesis representations zi1i2i3 with observation number zero being the gap filler observation. Equation (0.1)(b) requires that each of the tracks T1, . . . , TN1 be extended by one and only one hypothesis. Equation (0.1)(c) relates to each point or observation in the first observation set and requires that each such observation, except the gap filler, can only associate with one track but, because of the xe2x80x9cless thanxe2x80x9d condition, it might not associate with any track. Equation (0.1)(d) is like (0.1)(c) except that it is applicable to the second observation set. Equation (0.1)(e) requires that elements of the solution matrix z be limited to the zero and one values.
The only known method to solve Problem Formulation (0.1) exactly is a method called xe2x80x9cBranch and Bound.xe2x80x9d This method provides a systematic ordering of the potential solutions so that solutions with a same partial solution are accessible via a branch of a tree describing all possible solutions whereby the cost of unexamined solutions on a branch are incrementally developed as the cost for other solutions on the branch are determined. When the developing cost grows to exceed the previously known minimal cost (i.e., the bound) then enumeration of the tree branch terminates. Evaluation continues with a new branch. If evaluation of the cost of a particular branch completes, then that branch has lower cost than the previous bound so the new cost replaces the old bound. When all possible branches are evaluated or eliminated then the branch that had resulted in the last used bound is the solution. If we assume that N1=N2=N3=n and that branches typically evaluate to half there full length, then workload associated with xe2x80x9cbranch and boundxe2x80x9d is proportional to (n!|n/2!)2. This workload is unsuited to real time evaluation.
The Branch and Bound algorithm has been used in past research on the Traveling Salesman Problem. Messrs. Held and Karp showed that if the set of constraints was relaxed by a method of Lagrangian Multipliers (described in more detail below) then tight lower bounds could be developed in advance of enumerating any branch of the potential solution. By initiating the branch and bound algorithm with such a tight lower bound, significant performance improvements result in that branches will typically evaluate to less than half their full length.
Messrs. Frieze and Yadagar in dealing with a problem related to scheduling combinations of resources, as in job, worker and work site, showed that Problem Formulation (0.1) applied. They further described a solution method based upon an extension of the Lagrangian Relaxation previously mentioned. The two critical extensions provided by Messrs. Frieze and Yadagar were: (1) an iterative procedure that permitted the lower bound on the solution to be improved (by xe2x80x9chill climbingxe2x80x9d described below) and (2) the recognition that when the lower bound of the relaxed problem was maximized, then there existed a method to recover the solution of the non-relaxed problem in most cases using parameters determined at the maximum. The procedures attributed to Messrs. Frieze and Yadagar are only applicable to the three-dimensional problem posed by Problem Formulation (0.1) and where the cost matrix is fully populated. However, tracking multiple airborne objects usually requires solution of a much higher dimensional problem.
FIGS. 1 and 2 illustrate an example where xe2x80x9clook aheadxe2x80x9d data from the second image improved the assignment accuracy for the first image. Without the look ahead, and based only upon a simple nearest neighbor approach, the assignments in the first set would have been reversed. Problem Formulation (0.1) and the prior art only permit looking ahead one image. In the prior art it was known that the accuracy of assignments will improve if the process looks further ahead, however no practical method to optimally incorporate look ahead data existed. Many real radar tracking problems involve hundreds of tracks, thousands of observations per observation set and matrices with dimensions in the range of 3 to 25 including many images of look ahead.
It was also known that the data assignment problem may be simplified (without reducing the dimension of the assignment problem) by eliminating from consideration for each track those points which, after considering estimated limits of speed and turning ability of the objects, could not physically be part of the track. One such technique, denoted hereinafter the xe2x80x9ccone method,xe2x80x9d defines a cone as a continuation of each previously determined track with the apex of the cone at the end of the previously defined track. The length of the cone is based on the estimated maximum speed of the object and the size of the arc of the cone is based on the estimated maximum turning ability of the object. Thus, the cone defines a region outside of which no point could physically be part of the respective track. For any such points outside of the cones, an infinite number could be put in the cost matrix and a zero could be preassigned in the assignment matrix. It was known for the tracking problem that these elements will be very common in the cost and selection matrices (so these matrices are xe2x80x9csparsexe2x80x9d).
It was also known in the prior art that one or more tracks which are substantially separated geographically from the other tracks can be separated also in the assignment problem. This is done by examining the distances from each point to the various possible tracks. If the distances from one set of points are reasonably short only in relation to one track, then they are assigned to that track and not further considered with the remainder of the points. Similarly, if a larger group of points can only be assigned to a few tracks, then the group is considered a different assignment problem. Because the complexity of assignment problems increases dramatically with the number of possible tracks and the total number of points in each matrix, this partitioning of the group of points into a separate assignment problem and removal of these points from the matrices for the remaining points, substantially reduces the complexity of the overall assignment problem.
A previously known Multiple Hypothesis Testing (MHT) algorithm (see Chapter 10 of S. S. Blackman. Multiple-Target Tracking with Radar Applications. Artech House, Norwood, Mass., 1986.) related to formulation and scoring. The MHT procedure describes how to formulate the sparse set of all reasonable extension hypothesis (for FIG. 1 the set {h11 . . . h24}) and how to calculate a cost of the hypothesis {Ti, o1j, o2k} based upon the previously calculated cost for hypothesis (Ti, o1j}. The experience with the MHT algorithm, known in the prior art, is the basis for the assertion that look ahead through k sets of observations results in improved assignment of observations from the first set to the track.
In theory, the MHT procedure uses the extendable costing procedure to defer assignment decision until the accumulated evidence supporting the assignment becomes overwhelming. When it makes the assignment decision, it then eliminates all potential assignments invalidated by the decision in a process called xe2x80x9cpruning the tree.xe2x80x9d In practice, the MHT hypothesis tree is limited to a fixed number of generations and the overwhelming evidence rule is replaced by a most likely and feasible rule. This rule considers each track independently of others and is therefore a local decision rule.
A general object of the present invention is to provide an efficient and accurate process for assigning each point object in a region from multiple images to a proper track and then taking an action based upon the assignments.
A more specific object of the present invention is to provide a technique of the foregoing type which determines the solution of a k-dimensional assignment problem where xe2x80x9ckxe2x80x9d is greater than or equal to three.
The present invention relates to a method and apparatus for tracking objects. In particular, the present invention tracks movement or trajectories of objects by analyzing radiation reflected from the objects, the invention being especially useful for real-time tracking in noisy environments.
In providing such a tracking capability, a region containing the objects is repeatedly scanned to generate a multiplicity of sequential images or data observation sets of the region. One or more points (or equivalently observations), in each of the images or observation sets are detected wherein each such observation either corresponds to an actual location of an object or is an erroneous data point due to noise. Subsequently, for each observation detected, figures of merit or costs are determined for assigning the observation to each of a plurality of previously determined tracks. *Afterwards, a first optimization problem is specified which includes:
(a) a first objective function for relating the above mentioned costs to potential track extensions through the detected observations (or simply observations); and
(b) a first collection of constraint sets wherein each constraint set includes constraints to be satisfied by the observations in a particular scan to which the constraint set is related. In general, there is a constraint for each observation of the scan, wherein the constraint indicates the number of track extensions to which the observation may belong.
In particular, the first optimization problem is formulated, generated or defined as an M-dimensional assignment problem wherein there are M constraint sets in the first collection of constraint sets (i.e., there are M scans being examined) and the first objective function minimizes a total cost for assigning observations to various track extensions wherein terms are included in the cost, such that the terms have the figures of merit or costs for hypothesized combinations of assignments of the observations to the tracks. Subsequently, the formulated M-dimensional assignment problem is solved by reducing the complexity of the problem by generating one or more optimization problems each having a lower dimension and then solving each lower dimension optimization problem. That is, the M-dimensional assignment problem is solved by solving a plurality of optimization problems each having a lower number of constraint sets.
The reduction of the M-dimensional assignment problem to a lower dimensioned problem is accomplished by relaxing the constraints on the points of one or more scans thereby permitting these points to be assigned to more than one track extension. In relaxing the constraints, terms having penalty factors are added into the objective function thereby increasing the total cost of an assignment when one or more points are assigned to more than one track. Thus, the reduction in complexity by this relaxation process is iteratively repeated until a sufficiently low dimension is attained such that the lower dimensional problem may be solved directly by known techniques.
In one embodiment of the invention, each k-dimensional assignment problem (2 less than k less than M) is iteratively reduced to a (kxe2x88x921)-dimensional problem until a two-dimensional problem is specified or formulated. Subsequently, the two-dimensional problem formulated is solved directly and a xe2x80x9crecoveryxe2x80x9d technique is used to iteratively recover an optimal or near-optimal solution to each k-dimensional problem from a derived (kxe2x88x921)-dimensional problem k=2,3,4, . . . M.
In performing each recovery step (of obtaining a solution to a k-dimensional problem using a solution to a (kxe2x88x921)-dimensional problem), an auxiliary function is utilized. In particular, to recover an optimal or near-optimal solution to a k-dimensional problem, an auxiliary function, xcexa8kxe2x88x921 (k=4,5, . . . ,M), is specified and a region or domain is determined wherein this function is maximized, whereby values of the region determine the penalty factors of the (kxe2x88x921)-dimensional problem such that another two-dimensional problem can be formulated which determines a solution to the k-dimensional problem using the penalty factors of the (kxe2x88x921)-dimensional problem.
Each Auxiliary function xcexa8k, k=3, . . . ,Mxe2x88x921, is a function of both lower dimensional problems, penalty factors, and a solution at the dimension k at which the penalized cost function is solved directly (typically a two-dimensional problem). Further, in determining, for auxiliary function xcexa8k, a peak region, a gradient of the auxiliary function is determined, and utilized to identify the peak region. Thus, gradients are used for each of the approximation functions xcexa83,xcexa84, . . . ,xcexa8Mxe2x88x921 which are sequentially formulated and used in determining penalty factors until Mxe2x88x921 is used in determining the penalty factors for the (Mxe2x88x921)-dimensional problem. Subsequently, once the M-dimensional problem is solved (using a two-dimensional problem to go from an (Mxe2x88x921)-dimensional solution to an M-dimensional solution), one or more of the following actions are taken based on the track assignments: sending a warning to aircraft or a ground or sea facility, controlling air traffic, controlling anti-aircraft or anti-missile equipment, taking evasive action, working on one of the objects.
According to one feature of this first embodiment of the present invention, the following steps are also performed before the step of defining the auxiliary function. A preliminary auxiliary function is defined for each of the lower dimensional problems having a dimension equal or one greater than the dimension at which the penalized cost function is solved directly. The preliminary auxiliary function is a function of lower order penalty values and a solution at the dimension at which the penalized cost function was solved directly. In determining a gradient of the preliminary auxiliary function, step in the direction of the gradient to identify a peak region of the preliminary auxiliary function and determine penalty factors at the peak region. Iteratively repeat the defining, gradient determining, stepping and peak determining steps to define auxiliary functions at successively higher dimensions until the auxiliary function dimension (kxe2x88x921) is determined. In an alternative second embodiment of the present invention, instead of reducing the dimentionality of the M-dimensional assignment problem by a single dimension at a time, a plurality of dimensions are relaxed simultaneously. This new strategy has the advantage that when the M-dimensional problem is relaxed directly to a two-dimensional assignment problem, then all computations may be performed precisely without utilizing an auxiliary function such as xcexa8k as in the first embodiment. More particularly, the second embodiment solves the first optimization problem (i.e., the M-dimensional assignment problem) by specifying (i.e., creating, generating, formulating and/or defining) a second optimization problem. The second optimization problem includes a second objective function and a second collection of constraint sets wherein:
a) the second objective function is a combination of the first objective function and penalty factors or terms determined for incorporating (Mxe2x88x92m)-constraint sets of the first optimization problem into the second objective function;
b) the constraint sets of the second collection include only m of the constraint sets of the first collection of constraints,
wherein 2 less than m less than Mxe2x88x921. Note that, once the second optimization problem has been specified or formulated, an optimal or near-optimal solution is determined and that solution is used in specifying (i.e, creating, generating, formulating and/or defining) a third optimization problem of (Mxe2x88x92m) dimensions (or equivalently constraint sets). The third optimization problem is subsequently solved by decomposing it using the same procedure of this second embodiment as was used to decompose the first optimization problem above. Thus, a plurality of instantiations of the third optimization problem are specified, each successive instantiation having a lower number of dimensions, until an instance of the third optimization problem is a two-dimensional assignment problem which can be solved directly. Subsequently, whenever an instance of the third optimization problem is solved, the solution is used to recover a solution to the instance of the first optimization problem from which this instance of the third optimization was derived. Thus, an optimal or near-optimal solution to the original first optimization problem may be recovered through iteration of the above steps.
As mentioned above, the second embodiment of the present invention is especially advantageous when m=2, since in this case all computations may be performed precisely and without utilizing auxiliary functions.
Other features and benefits of the present invention will become apparent from the detailed description with the accompanying drawings contained hereinafter.